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Abstract. Computer simulations of simple model systems for liquid crystals are briefly 
reviewed, with special emphasis on systems of ellipsoids. First, we give an overview 
over some of the most commonly studied systems (ellipsoids, Gay-Berne particles, 
spherocylinders) . Then we discuss the structure of the nematic phase in the bulk and 
at interfaces. 



1.1 Introduction 

Randomly distributed hard spheres in three dimensions form two types of struc- 
tures, depending on their density: Fluid and crystalline. For randomly dis- 
tributed anisotropic particles, the situation can be different: Several phases may 
exist at intermediate densities between the fluid and the solid state. These phases 
are called "mesophases" . For example, one often observes a nematic phase where 
the particles are oriented in one common preferred direction, but have no crys- 
talline translational order. Other common phases are the smectic phases, in 
which the particles are arranged in layers of two dimensional fluids. Some of 
these structures are sketched in Figure 1. Since the mesophases are neither crys- 
talline nor truly liquid, they are commonly referred to as "liquid crystal phases" . 



Isotropic I 



3 



Nematic N 



Smectic A 
(Smectic B 



'ci:;> 

Sm A 

) Sm BHex 



Smectic C Sm c 
(Smectic F) Sm f 



Figure 1: 

Liquid crystalline phases. In the Smectic Bjjg-j^^ and Smectic F phase, the structure 
within the layers is not entirely fluid, but has a type of order called hexatic. See Ref. 
plH for explanation. 



Experimentally, such phases are observed in various systems of anisotropic parti- 
cles: Small organic molecules, stiff polymers, self-associated wormlike micelles, or 
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even whole virusses (one famous example being the tobacco mosaic virus). The 
phase transitions are triggered in some cases by the temperature (thermotropic 
liquid crystals), and in other cases by the concentration of the anisotropic parti- 
cles (lyotropic liquid crystals). In reality, particles are of course not distributed 
"randomly", but according to a Boltzmann distribution P oc e~^/^^'^ which 
accounts for the energy E of their mutual interactions (T is the temperature 
and Ub the Boltzmann constant). Nevertheless, the mechanisms which drive the 
phase transitions in lyotropic liquid crystals are essentially the same as those in 
systems of hard anisotropic particles. Therefore, the latter are often employed to 
model phase transitions in liquid crystals. More generally, systems of anisotropic 
particles with simplified interactions have proven very useful to study generic 
properties of liquid crystals • 

In the present contribution, we describe some of the computer simulation work 
that has been done in this direction. We do not intend to give a complete account 
of this active and rapidly progressing field of research. The limited space here 
permits only a rather crude introduction and the discussion of a few examples. 
The reader who is interested in more exhaustive overviews is referred to, e.g. , 
the set of review articles in Ref. [§, or to Refs. 

The paper is organized as follows: First, we review some of the most commonly 
studied idealized models for particle based simulations. We will then focus on 
the nematic phase and discuss some of its bulk properties. Finally, we address 
specific issues of intcrfacial properties in nematic liquid crystals. 



1.2 Model systems 



1.2.1 Ellipsoids 



Perhaps the most obvious anisotropic generalization of hard spheres are hard 
ellipsoids of revolution with one symmetry axis of length L and transverse thick- 
ness D. The phase diagram in three dimensions has been established from com- 



puter simulations by Frenkel, Allen and coworkers [yjlC|,|ll|,|l2 13 IJ]. It is shown 
for a range of elongations k — L/D in Figure 2. If the particles arc sufficiently 
anisotropic, a nematic phase (N) intrudes between the isotropic (I) and the 
crystalline (X,PC) phase. No further liquid crystalline phases are present. In 
particular, smectic phases do not exist. This makes hard ellipsoids particularly 
suited to study the properties of nematic liquid crystals. 

From a technical point of view, the study of ellipsoids is complicated by the 
fact that the contact distance cr(ui, Uj, ) between two particles i and j with 
orientations Uj and in the center-center direction Yij {\vLij\ = |fij| — 1)) 
cannot be given in analytically closed form. Criteria to decide whether or not 
two particles overlap have been derived by Vieillard-Baron [ p^ , and by Perram 
and Wertheim Technical details on the implementation can be found in 

the review article Ref. |Tl|]). 
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Figure 2: 

Piiase diagram of liard ellipsoids as a function of elongation k — L/D and reduced 
density p/ pep {pep being the density at hep packing). Phases are: I (isotropic phase), 
N (nematic phase), X (orientationally ordered crystalline phase), PC (orientationally 
disordered crystalline phase), (a) Full range of elongations: From Ref. with data 
from Refs. ||jl0|jl|Jl|. (b) Detail from Ref. Points denote simulation data, lines 
predictions from different theories. 



Both the simulation and the data analysis are simplified considerably if an ex- 
plicit formula for <7{ui , Uj , f ij ) is available. For this reason, Berne and Pechukas 
have proposed an approximate expression 



cr(Ui 



Uj,r,j) = (To 
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1 - XU» • Uj 

(1-1) 

with the anisotropy parameter x = {i^^ ~ + 1) This function gives the 

exact contact distance if the particles i and j are co-linear (parallel to each 

< \/2 



other), and overestimates it by a factor of at most \/2/ v 1 + 2k/ (k^ -|- 1) 
otherwise. A few simulation studies have considered systems of "hard" parti- 
cles with contact distances given by Eqn. (1.1) ("hard Gaussian overlap parti- 
cles") |20|^ , p^^3y24|25|| . Much more commonly, however, the contact function 
( [l.l| ) is used in conjunction with Lcnnard- Jones type interaction potentials, the 
Gay-Berne potential . This class of models will be discussed in the next sec- 
tion. The present authors favor yet another potential, which is purely repulsive 
but soft. 



V,.. 



0, 



otherwise 



(1.2) 



where Ra is a reduced distance, 



Ri 



{vij - cr(ui,Uj,f,j) +a-o)/cro 



(1.3) 



As far as we have seen so far, systems of such "soft ellipsoids" (at temperature 
ksT ^ e) have almost the same structure as systems of hard ellipsoids. 
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Figure 3: 

Phase diagram of Gay-Berne ellipsoids with fi — 2, u = 1, k — 4, and k' = 5 as a 
function of temperature T in units of eo/ks and number density p . In addition to the 
phases observed for hard ellipsoids (Figure 2) , one also finds smectic phases Sa and Sb ■ 
The simulations did not allow to decide whether the smectic B phase is hexatic or solid, 
i.e., crystalline. In the smectic A phase, however, the smectic layers are undoubtedly 
fluid. From Ref. Q 



1.2.2 Gay-Berne particles 



Despite their appealing simplicity, hard or soft repulsive ellipsoids are not in 
every respect best suited to study generic liquid crystal properties. For example, 
such models do not exhibit smectic phases. Furthermore, one cannot study the 
effect of attractive interactions, the interplay between liquid-crystalline order 
and liquid-vapor phase separation etc. Therefore, Gay and Berne have in- 
troduced a class of ellipsoidal pair potentials with attractive interactions, which 
has become one of the standard liquid crystal potentials. 



The functional form of these potentials is 



(1.4) 



where the reduced distance Rij defined is as above (eqn. (1.3) with <j taken from 
(pTl|)), and the energy functions e are given by 



e'(u,,Uj,f,j) = 1 - — 



e"(Uj,Uj 



1 + x'Uj • Uj 

[l-x^(u,.u,f]-^/^ 



(U^ ■ f,j- - VLj ■ TjjY 

1 - x'Ui ■ Uj 



(1.5) 
(1.6) 
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The new additional anisotropy parameter x' accounts for the fact that the at- 
tractive energy should favor a side-by-side alignment of particles compared to an 
end-to-end alignment. It is related to the corresponding ratio of minimum ener- 
gies (well depths) k! ^ e^ido-^ido/eond-cnd via x' = W^^^ - + 1). Gay and 
Berne originally suggested to choose the exponents /i = 2 and 1^ = 1. The phase 
behavior for this set of exponents has been studied in detail p^ , |28| , p9| , |30y3 1|J3^ , 
for various choices of k and k' [ p0pl| , |3^ . An example is shown in Figure 3. One 
observes a liquid-vapor coexistence region at low temperature. Furthermore, the 
attractive interaction stabilizes a smectic phase (Sm A), which consists of stacked 
layers of two dimensional fluids. 

The Gay-Berne model is widely used in liquid crystal simulations |^ . Due to the 
option of varying not only k and k', but also the exponents ^ and ly, it is very 
versatile jS^^Q. Luckhurst et al have demonstrated that it can even be adjusted 
to serve as a good model for a real thermotropic liquid crystal 



35, 3f 



1.2.3 Other models 

Ellipsoids are not the only particles that have been used to model liquid crys- 
tals on an idealized level. An almost equal amount of work has been devoted to 
systems of spherocylinders: Cylinders of length L and diameter D capped with 
hemispheres at both ends. They have the advantage that the contact distance 
a between two particles can be evaluated exactly pTf . Furthermore, they re- 
semble quite closely certain real colloidal liquid crystal substances, e. g., rodlike 
virusses 



The phase behavior of hard spherocylinders has been studied by Frenkel and 
others [ pp^ , ^ , ^5p^ , p7t . The complete phase diagram for all elongations L /D 
as computed by Bolhuis and Frenkel |Q is shown in Figure 4. It differs strikingly 
from the corresponding phase diagram for hard ellipsoids in that sufficiently 
elongated spherocylinders form smectic phases. This came as a surprise when it 
was first discovered in simulatio ns |4^ . Meanwhile, it has been reproduced by 
density functional theories [ p| , p|j5C| , |5l| , p2[ . 

Bolhuis et al have studied the influence of attractive interactions on the phase 
behavior of spherocylinders. Like in systems of Gay-Berne particles (Figure 3), a 
liquid- vapor coexistence region emerges at low temperatures (see also Ref. js^). 
The influence of polydispersity - i. e., fluctuating rod length - on the phase 
diagram has been investigated by Bates et al |5^. Perhaps not surprisingly, the 
smectic phase is destabilized in favor of the nematic phase. 

We shall only briefly mention some other particle based model systems for liquid 
crystals: In a number of studies, liquid crystalline molecules have been modeled 
by stiff chain molecules |5^ , |57| , p8| , p9|j60| , |6l| , |62| , ^ . with aspect ratio k < 1 and 
spheres cut at both sides have been used to study discotic liquid crystals, i. e., 
fluids of disklike particles ^4 65 6^ , |67| , ^j69| . Such fluids exhibit new, different 



types of liquid crystalline order. For example, the smectic phases are replaced 
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by columnar phases: The discs assemble into hexagonal arrays of columns, but 
remain fluid in one dimension within the columns. 
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Figure 4: 

Phase diagram of hard spherocylinders as a function of reduced density p* = p/ pep 
and elongation L/D (left) or inverse elongation D/L (right). Phases are: I (isotropic 
phase), N (nematic phase), Sm (smectic A phase), AAA, S=ABC, (two types of solid 
phases with orientational order), and P (orientationnally disordered solid phase). From 
Ref. @. The transition between the nematic and the smectic phase is first order at 
all elongations |47|]. 



1.3 Properties of the nematic phase 

A number of important properties of the nematic phase can already be deduced 
from very general considerations. 

First, geometric arguments suggest that the transition between the isotropic and 
the nematic phase should be discontinuous or first order, i.e., there should be a 
density regime in which both phases coexist. This is a consequence of "Landau's 
rule" , which relates the character of the phase transition to the symmetry 

groups of the two phases. Landau's rule is not rigorous, it does not exclude the 
possibility that the width of the coexistence region may happen to shrink to 
zero in certain points of the phase diagram. Moreover, order fluctuations may 
affect the nature of the phase transition. So far, however, the actually observed 
isotropic-nematic transitions have been discontinuous, albeit with rather narrow 
coexistence regions. 

Second, the free energy cost of spatial modulations of the average particle ori- 
entation must vanish if the wavelength of the modulations tends to infinity. 
This is an example of Goldstone's theorem [ [ziHt^ : If the order parameter in 
an ordered phase breaks a continuous symmetry, there exist soft massless fluc- 
tuation modes |7^. Modulations with finite wavelength are subject to elastic 
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restoring forces |73[ , |74| |. For symmetry reasons, these depend on only three ma- 
terial constants, the Frank elastic constants Kn, K22, K33 [^5|j7^,^. Figure 5 
illustrates the corresponding fundamental distortions, the "splay", "twist" and 
"bend" mode. 




Splay Twist Bend 



Figure 5: 

Elastic modes in nematic liquid crystals. 



We turn to give a more quantitative description. The nematic order is usually 
characterized by the 3x3 dimensional nematic order tensor 

1^3 1 

1=1 

where the sum runs over all particles i, 1 denotes the unit matrix and (g) the 
dyadic product. The largest eigenvalue of Q is the nematic order parameter 
per particle, and the corresponding eigenvector is the director n of the nematic 
liquid. It gives the direction of preferential alignment of the molecules. Long 
wavelength spatial modulations n(r) of the director are penalized by the elastic 
free energy Q 

.F{n(r)} = i y drji^n [V • n]^ + K22[n ■ (V x n)]^ + Kssin x (V x n)]^}. (1.8) 

On macroscopic length scales, the structure of liquid crystals is controlled almost 
exclusively by the elastic energy and thus by the Frank elastic constants ifn, 
K22, and K33. 

The quantitative characterization of the local, microscopic structure is much 
more involved. The most intensely studied quantities are usually the pair distri- 
butions or the pair correlation functions. On principle, a hierarchy of infinitely 
many A'^-particle distribution functions would be required to fully characterize a 
fluid. In systems with only pairwise interactions, however, the pair distributions 
stand out because they can be used to calculate several macroscopic quanti- 
ties such as the pressure and the compressibility. Moreover, pair correlations are 
often accessible to experiments in real systems. 
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As before, we characterize the orientation of a particle by a unit vector u. The 
one-particle distribution p^^^ (u, r) gives the probability density of finding a par- 
ticle of orientation u at the position r. Similarly, the two-particle distribution 
p^'^\ui,ri,Uj,rj) is the probability density of finding simultaneously one parti- 
cle of orientation at the position r^, and another with orientation Uj at the 
position Tj. In a homogeneous nematic phase, p^^^ does not depend on the posi- 
tion r, and p^^^ depends only on the relative position r^j = — r^. Furthermore, 
particles at infinite distances are uncorrelated, hence 

p(2)(u„u„r)''^/i)(ui)p«(u2). (1.9) 

The correlations between particles at finite distance are characterized by the 
total correlation function ||7^ , [7^ 

h{Ul,U2,r)^-rj,-—y——--l, 1.10 

It subsumes the total effect of a particle i on a particle j. Due to the elastic 
interactions mentioned above, it decays only slowly at long distances with the 
algebraic power law l/r. 

Of course, other definitions of correlation functions are possible. As an alternative 
to the total correlation function, one could consider the orientation correlation 
function defined by Penttinen and Stoyan . It describes the pure orientation 
correlations between particles at a given distance, and should display the same 
asymptotic l/r power law in the nematic fluid as h. 

The long range correlations are mediated by the whole bulk of the surrounding 
nematic fluid. It seems desirable to separate these "indirect" effects from a more 
local "direct" effect of two particles on each other. This can be done by con- 
sidering the "direct correlation function" c(u,i, Uj, r^), which is defined through 
the Ornstein-Zernike equation ||78| , [79| 

h{u,,Uj,rij) = c{u,,Uj,rij) + J c{ui,Uk,r,k) p^^H^k) h{uk,Uj,rkj)dukdrk■ 
il.n) 

As it turns out, the direct correlation function is indeed short range even in the 
nematic phase. 

We illustrate this with recent computer simulation results. We have performed 
computer simulations of 1000-8000 ellipsoids with elongation k = 3 at the num- 



ber density p = 0.3, interacting with the potential (1.2) at temperature T = 0.5. 
Technical details can be found in Refs. [plU 



The calculation of direct correlation functions from computer simulations is 
computationally demanding and involves a complicated data analysis. We were 
the first to determine the direct correlation function from computer simulations 
of a nematic fluid without any approximations. Direct correlation functions in 
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isotropic fluids have been calculated earlier by Allen et al [ p5pq |, and approx- 
imate data for ncmatic fluids have been derived from simulations by Stelzer et 
al[|g|8|. 

Figure 6 shows the orientational average of the total correlation function and the 
direct correlation function. The average is carried out over the orientations u^, Uj, 
and Tij. Figure 6 demonstrates that long range algebraic correlations between 
particles are not apparent in this averaged curve. The elastic interactions affect 
only orientations, not total densities. Therefore, both h and c are short range. 
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Figure 6: 

Total correlation function h (left) and direct correlation function c (right) vs. molecular 
distance r, averaged over all orientations of Ui,Uj, and r^j, for different system sizes 
N. 



In order to assess the effect of the elasticity, we need to study orientational de- 
pendent correlations. To this end, we expand all correlation functions in spherical 
harmonics Yim{u). 

(r) Yl^miiUi) ^2^2 (Uj) Yijni^), (1-12) 



where F stands for p^'^\ h, or c. The z axis is chosen in the direction of the direc- 
tor. The orientation dependence of the correlations is reflected by the coefficients 
with nonzero indices. Figure 7 (left) shows a coefficient of the total correlation 
function with a particularly pronounced long range tail. As shown in Figure 7 
(right), the tail completely disappears in the direct correlation function. 

The direct correlation function is a central quantity in theories of liquid mat- 



ter 78 1 . It enters density functional theories which predict the structure of in- 
homogeneous fluids, e. g., fluids at interfaces. It allows to calculate a number 
of material constants of fluids. In nematic fluids, for example, it can be used 
to calculate the Frank elastic constants according to a set of equations flrst de- 
rived by Poniewierski and Stecki [p9lM . We have applied these equations to our 
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data and obtained a set of values for Kn, K22, and which is in agreement 
with results from a different, more direct method pT| , p2[ . Thus the direct corre- 
lation function establishes a bridge between the microscopic structure and the 
parameters which characterize the macroscopic structure. 
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Figure 7: 

Expansion coefficient with h = 2,l2 — 2, mi = l,m2 = —1,/ ~ 2, and m = of 
the total correlation function h (left), and the direct correlation function (right) vs. 
molecular distance r for different system sizes A'^. Left Figure: Dashed line indicates an 
extrapolation with a 1/r power law behavior. Inset shows the same data vs. 1/r. 



1.4 Interfacial properties 



We turn to discuss briefly the interfacial properties in nematic liquid crystals. 
These are of great interest in the liquid crystal display technology There- 
fore, a growing number of simulatio ns are devoted to the study of mo del nemat- 
ics at surfaces jo^ Hp5p^ ,|97|,p8|,pl |l00| , [l0l| |l 02yi03|,[r0^,|l05yi 06|]l 07|| . and inter- 



faces [|l08|]l09l|TT0| , |lll| , |ll2j or in thin films |^4|,|2 5|, |l 1 3l|l 1 4| , |l 1 5| , |l 1 GUi T7| , [iT8 | 



The microscopic structure at surfaces can be quite complex [^3|j9j,|95|,^ . From 
a macroscopic point of view, the presence of surfaces or interfaces introduces 
essentially three new effective parameters: the interfacial tension, the anchoring 
angle, and the anchoring energy. The anchoring angle is the angle with respect 
to the surface which the director preferably assumes close to the surface. The 
anchoring energy is related to the force needed to twist the director out of the 
anchoring angle. 

The anchoring angle can be calculated from simulations in a straightforward way. 
Obtaining the anchoring energy is more difficult and time consuming. Allen and 
coworkers have devised difllerent methods and applied them to systems of ellip- 
soids at hard, planar walls | IOC, 107 1. The interfacial tension can be computed 



from the anisotropy of the pressure tensor close to the surface 111 119| , 



the case of fiuid-fiuid interfaces from their undulations (capillary waves) |112| 
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Recently, a number studies have been devoted to the interface between the ne- 
matic and the isotropic phase in Gay-Berne fluids \10i] and fluids of spherocylin- 
ders [|ll§ or ellipsoids ]l09| , pTl|pT2| . They reveal among other a rather intriguing 
capillary wave spectrum, ft reflects a complex interplay between the bare surface 
tension and the elastic interactions in the ncmatic phase |112|. 



1.5 Conclusions 



In materials science, systems of ellipsoids serve as model systems which exhibit 
generic properties of nematic liquid crystals. They are studied for the general 
insight they can give into the physics of liquid crystals, and also because sim- 
ulations of ellipsoids are computationally much cheaper than simulations of a 
comparable number of more realistic molecules. However, fluids of ellipsoids are 
also fascinating from a much more general point of view: They are simple systems 
where purely entropic effects generate an amazing variety of different structures, 
and lead to unusual long-range static and dynamic properties. Computer simu- 
lations are often the only way of investigating some of that behavior. 
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